
require(Rcpp)
setwd(wkdir)
output_version <- '2024-07-19'
datdir <- file.path(wkdir,'local')

rasterOptions(tmpdir = "C:/temp")

fishnet_lake_chad_ntl_adm0_x_d01_shp <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01')
fishnet_adm0_only <- read_sf(dirname(fishnet_lake_chad_ntl_adm0_x_d01_shp),basename(fishnet_lake_chad_ntl_adm0_x_d01_shp)) %>%
  dplyr::select('OBJECTID','GID_0')
st_crs(fishnet_adm0_only) <- "EPSG:4326"

# raster of grid dd 0.1
r_lakechad <- terra::rast(extent = c(0,25,0,25),res=0.1, crs="EPSG:4326")
r_fishnet_adm0_only <- rasterize(fishnet_adm0_only,r_lakechad,field="OBJECTID")

##
## class 10 = cropland
## class 20 = cropland / irrigated

lulc08_fn <- file.path(datdir,'gis_data/007_environment/lulc/lulcre_factor2008.tif')

if(!file.exists(lulc08_fn)){
  lulc9215 <- stack(paste0(lulc_dir,'/ESACCI-LC-L4-LCCS-Map-300m-P1Y-1992_2015-v2.0.7.tif'))
  names(lulc9215)<- paste0("lulc",1992:2015)
  
  lulc08 <- lulc9215[["lulc2008"]]
  
  # manually set outside for 5x5
  ntl.extent <- c(0,25,0,25)
  
  # crop extent to lake chad area
  lulc = crop(lulc08,ntl.extent)
  
  # reclass matrix
  # 0 = no data
  # 1 = cropland
  # 2 = cropland mosaic
  # 3 = grassland (130)
  # 4 = urban
  # 5 = bare 
  # 6 = other
  # 7 = water
  # "cropland", "mosaic", "other", "urban" #
  m <- c(0,0,1,10,6,11,6,12,1,20,2,30,6,40,6,50,6,60,6,61,6,62,6,70,6,71,6,72,6,80,6,81,6,
         82,6,90,6,100,6,110,6,120,6,121,6,122,3,130,6,140,6,150,6,151,6,152,6,153,6,160,6,170,
         6,180,4,190,5,200,5,201,5,202,7,210,5,220)
  rclmat <- matrix(m, ncol=2, byrow=TRUE)
  # wrong order fix later #
  rclmat <- cbind(rclmat[,1],rclmat[,2])
  
  # set up the cluster object for parallel computing
  beginCluster(12)
  system.time(rclulc <- clusterR(lulc,reclassify,args=list(rclmat)))
  endCluster()
  
  names(rclulc)<-paste0("lulcre",2008)
  
  is.factor(rclulc)
  f_rclulc <- as.factor(rclulc)
  labels0 <- c("cropland","cropland_mosaic","grassland",
               "urban","bare","other","water")
  # Specify your categories; you can customize this based on your needs
  rat <- levels(f_rclulc)[[1]]
  rat$lulc <- labels0
  levels(f_rclulc) <- rat
  #f_rclulc_fn <- paste0("lulcre_factor",2008)
  writeRaster(f_rclulc,
              filename=lulc08_fn, 
              bylayer=TRUE, driver='GTiff',overwrite=TRUE)
  
} else {
  rclulc <- rast(lulc08_fn)
  f_rclulc <- as.factor(rclulc)
  rat <- levels(f_rclulc)[[1]]
  labels0 <- c("cropland","cropland_mosaic","grassland",
               "urban","bare","other","water")
  # Specify your categories; you can customize this based on your needs
  rat <- levels(f_rclulc)[[1]]
  rat$lulc <- labels0
  levels(f_rclulc) <- rat
  
}



sum_cover <- function(x){
  list(x %>%
         group_by(value) %>%
         summarize(total_area = sum(coverage_area)) %>%
         mutate(sh_lulc = total_area/sum(total_area)))
  
}

##
## RADIUS BUFFER 0 (grid cell only) 50KM (FROM POINT), 100KM (FROM POINT)
##

## memory issues
n_objectid_list=unique(fishnet_adm0_only$OBJECTID)
total_chunks <- split(n_objectid_list, ceiling(seq_along(n_objectid_list)/1000))

require(doParallel)

for(chunk_index in 1:length(total_chunks)) {
  print(chunk_index)
  
  lulc_chunk_out_fn <- file.path(wkdir,'proc_data','lulc_chunks',sprintf('LULC08_share_%s_of_%s_%s.csv',chunk_index,length(total_chunks),output_version))
  xfun::dir_create(dirname(lulc_chunk_out_fn))
  if(!file.exists(lulc_chunk_out_fn)){
    n_objectid_list_chunk <- as.numeric(unlist(total_chunks[chunk_index]))
    
    n_cores <- min(16,length(n_objectid_list_chunk))
    
    cl <- makeCluster(n_cores)
    registerDoParallel(cl)
    
    results <- foreach(objectid_j=n_objectid_list_chunk,
                       .export=c(ls()),
                       .packages = c("dplyr", "sf","exactextractr",
                                     "spatialEco", "terra")) %dopar% {
                                       # Modulus operation
                                       if(objectid_j %% 100==0) {
                                         # Print on the screen some message
                                         cat(paste0("iteration: ", objectid_j, "\n"))
                                       }
                                       
                                       #udf
                                       sum_cover <- function(x){
                                         list(x %>%
                                                group_by(value) %>%
                                                summarize(total_area = sum(coverage_area)) %>%
                                                mutate(sh_lulc = total_area/sum(total_area)))
                                         
                                       }
                                       
                                       fishnet_lake_chad_ntl_adm0_x_d01_shp_fn <- paste0(wkdir,'/local/gis_data/003_boundaries/fishnet/fishnet_lake_chad_ntl_adm0_x_d01.shp')
                                       fishnet_adm0_only_j = read_sf(fishnet_lake_chad_ntl_adm0_x_d01_shp_fn, query = sprintf('SELECT * from fishnet_lake_chad_ntl_adm0_x_d01 WHERE OBJECTID = %s',objectid_j))
                                       st_crs(fishnet_adm0_only_j) <- "EPSG:4326"
                                       
                                       if(dim(fishnet_adm0_only_j)[1]>0){
                                         
                                         fishnet_adm0_only_buf100km_bbox <- round(st_bbox(spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=150000)),1)
                                         
                                         ## crop raster to 100km buffer
                                         f_rclulc <-  terra::rast(lulc08_fn)
                                         f_rclulc_crop <- terra::crop(f_rclulc,fishnet_adm0_only_buf100km_bbox,snap="out")
                                         
                                         ## 0km
                                         x <- exact_extract(f_rclulc_crop, fishnet_adm0_only_j, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
                                         
                                         #merge the list elements into a df
                                         lulc_share0 <- bind_rows(x, .id = "objectid") %>%
                                           mutate(buffer_km = 0)
                                         
                                         # clean-up
                                         rm(x)
                                         
                                         # projected 50km
                                         fishnet_adm0_only_buf50km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=50000)
                                         x <- exact_extract(f_rclulc_crop, fishnet_adm0_only_buf50km, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
                                         
                                         #merge the list elements into a df
                                         lulc_share50 <- bind_rows(x, .id = "objectid") %>%
                                           mutate(buffer_km = 50)
                                         
                                         # projected 100km
                                         fishnet_adm0_only_buf100km <- spatialEco::geo.buffer(x=st_centroid(fishnet_adm0_only_j), r=100000)
                                         x <- exact_extract(f_rclulc_crop, fishnet_adm0_only_buf100km, coverage_area = TRUE, summarize_df = TRUE, fun = sum_cover)
                                         
                                         #merge the list elements into a df
                                         lulc_share100 <- bind_rows(x, .id = "objectid") %>%
                                           mutate(buffer_km = 100)
                                         
                                         out <- rbind(lulc_share0,lulc_share50,lulc_share100)
                                         
                                         # objectid id by iteration
                                         out <- out %>%
                                           mutate(objectid = objectid_j) %>%
                                           dplyr::select(-total_area)
                                         
                                         # clean-up
                                         rm(fishnet_adm0_only_buf100km_bbox,fishnet_adm0_only_buf50km,fishnet_adm0_only_buf100km)
                                         rm(x,lulc_share0,lulc_share50,lulc_share100)
                                         
                                         return(out)
                                       }
                                     }
    
    parallel::stopCluster(cl)
    
    system.time(lulc_dt <- data.table::rbindlist(results))
    lulc_df <- lulc_dt %>% as.data.frame()
    names(lulc_df)
    dim(lulc_df)
    
    data.table::fwrite(lulc_df,lulc_chunk_out_fn)
    
    # clean-up
    rm(lulc_df)
    rm(lulc_dt)
    
  }
  
  
}

# merge all chuncks
lulc_df_files <- sort(list.files(dirname(lulc_chunk_out_fn),full.names=T))
tables <- lapply(lulc_df_files, function(z) read.csv(z))
lulc_df_all <- data.table::rbindlist(tables)

#
rat <- rat %>%
  dplyr::select(ID,lulc)
lulc_out <- merge(lulc_df_all,rat,by.x='value',by.y='ID',all.x=T) %>%
  dplyr::select(objectid,sh_lulc,lulc,buffer_km) %>%
  arrange(objectid,lulc,buffer_km)

# check max obs. X objectid X (3) 
table(lulc_out$objectid)

lulc_out_fn <- file.path(wkdir,'proc_data',sprintf('LULC08_share%s.dta',output_version))

# Adding the label value
attr(lulc_out$objectid, "label") <- "objectid of grid"
attr(lulc_out$sh_lulc, "label") <- 'Share of area covered by LULC : ESA'
attr(lulc_out$lulc, "label") <- 'LULC category : ESA'
attr(lulc_out$buffer_km, "label") <- 'Buffer'
attr(lulc_out, "label") <- sprintf("bh_spillover_wd_prep_landcover_reclass2008_only_vector.R last run on %s",Sys.Date())
head(lulc_out)
summary(lulc_out)

## Save the data
haven::write_dta(lulc_out, lulc_out_fn)
